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AN EXPERIMENTAL COMPARISON OF SEVERAL APPROACHES 
TO THE LINEAR LEAST SQUARES PROBLEM 


Eugene J. Lefferts 
Daniel P. Muhonen 


ABSTRACT 


Several computational algorithms for obtaining least square solutions to a 
system of equations are made. The a’gorithms include both conventional and 
pseudo inversion methods. Comparisons of running time, computer storage and 
accuracy of recovery are made for each algorithm in both single and double 
precision. 
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AN EXPERIMENTAL COMPARISON OF SEVERAL APPROACHES 
TO THE LINEAR LEAST SQUARES PROBLEM 


SUMMARY 


Experiments were made using five procedures to solve the linear least squares problem. 
These were: 

(1) The Gauss-Jordan algorithm on the normal equations 

(2) An Andree algorithm for the Penrose pseudoinverse on the normal equations 

(3) A Gram-Schmidt orthogonalization pseudoinversion scheme on the normal 
equations 

(4) The Gram-Schmidt procedure applied to the original rectangular data matrix 

(5) Householder's algorithm for direct triangulation of the original data matrix 

The experiments involved the recovery of polynomial coefficients in both double and 
single precision. Computer running time and storage requirements are presented. 

The following conclusions and recommendations can be made as a result of these 
experiments : 

• Procedures (4) and (5) above, which avoid the formation of the normal 
equations, give considerably better results than do the other schemes. 

• Procedure (4) should be used in double precision when there are no storage 
or timing restrictions. It will handle situations of reduced rank. 

• Procedure (5) should be used in double precision when timing but not storage 
restrictions exist. It does not handle situations of reduced rank. 

• Procedures (1), (2), and (3) give better results in double precision than do 
procedures (4) and (5) in single precision. 

• Procedure (1) is most efficient in running time and storage requirements, but 
it should be avoided whenever possible because it cannot handle situations of reduced 
computational rank. 

• Procedure (2) gives somewhat better results than procedure (3), but it requires 
much more storage. 
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AN EXPERIMENTAL COMPARISON OF SEVERAL APPROACHES 
TO THE LINEAR LEAST SQUARES PROBLEM 


INTRODUCTION 

In solving the linear lease squares problem, complications arise due to the 
introduction of errors in the computation. The magnitudes of these errors vary 
considerably with the procedures used to obtain a least squares solution. One dif- 
ficulty that frequently arises is that, because of accumulated truncation error in the 
computation, only a few decimal figures in the original data matrix are of significance 
in the least squares process. Hence, if only the computationally significant portions 
of the matrix elements are considered, a matrix of maximal rank may become one 
of reduced rank. Unless the computational errors are accounted for in some way, a 
meaningless solution can result. 

Five approaches to the least squares problem were investigated experimentally. 
Using the normal equations, comparisons were made between the Gauss-Jordan in- 
version algorithm and two pseudoinversion schemes: an Andree algorithm and a 
Gram-Schmidt orthogonalization procedure. The advantage of implementing a pseudo- 
inverse is that the minimal norm solution can be found when the matrix is not of 
maximal rank. Both of these pseudoinverse algorithms use a-priori truncation error 
limits to determine the effective rank of the matrix. In addition, experiments were 
made using two approaches which arrive at a solution directly from the rectangular 
data matrix without forming the normal equations. The Gram-Schmidt pseudo - 
inversion scheme is compared here with the Householder algorithm which forms the 
Cholesky decomposition of A™ A directly from A. The effect of truncation error 
resulting from the formation of the normal equations is quite evident from these 
experiments. Computer running time and storage information are also presented in 
order that a more complete evaluation can be made of the advantages and disadvantages 
of all procedures. 


PROBLEM STATEMENT 

This presentation will be concerned with the linear system 

AX = B , 


where A is a known m x n matrix, B is a known m x 1 vector, _and X is an unknown 
solution vector n x 1. "*he system has a solution if and only if B is in the column 
space of A. If m > n, as is the case in many physical situations, this condition will 
generally not be satisfied. However, it is still possible to find absolution that is 
'best*' in the least squares sense; that is, one can find a vector X Q such that 

|| AX q - B || = min. . 

where 1 1 1 1 represents the vector Euclidean norm. If A is of maximal column rank 
(i. e. , rank (A) = n), then the solution X Q is unique, and 


1 


-* T — 1 T — • 

X 0 = (A a A) 1 A* b . 


This matrix equation represents what is commonly known as the normal equations. 

The solution X n in the above equation can be found by actually evaluating the 
u t 

square, nonsingular matrix A A and finding its inverse by some standard algorithm. 
There are computational disadvantages in forming the normal equations, however, 
since information may be lost by truncation in evaluating the inner products which 
form the elements of A T A (ref. 1). There is an alternative approach developed by 

Householder which can be used to find thf' solution vector X n without requiring the 

T ^ 

formation of A A (see Appendix C). In solving a system of n linear equations in n 

unknowns, a common procedure is to reduce the matrix of coefficients to upper tri- 
angular form by a series of elementary row operations (ref. 2)-, If these transfor- 
mations are performed simultaneously on the constant vector A 1 B, the solution 
vector can then be found by back substitution. In our problem this can be represented 
by the following equation: 

T -* T — 

P A AX 0 = P A B , 


where P represents the product of the elementary transformation matrices necessary 
to reduce A A to upper triangular form; that is, PA a A is an upper triangular matrix. 
Householder’s scheme provides for the evaluation of PA T A and PA^B directly from 
A and B without requiring the formation of the inner products which can result in 
serious truncation error. 


The problem still remains that the matrix A may not be of maximal rank, in 
which case the normal equations will not have a unique solution. From a computa- 
tional point of view, this is a very real problem in many physical situations. 

Certainly if the level of noise is sufficiently high in the data represented by two or 
more columns in the matrix A, very few figures of significance exist in this data. 
What significance there is can easily be lost in the least squares solution processes. 

It may also be true that a near-linear dependence does exist in the observables 
represented by two or more columns. Unless some means is provided to allow for 
these errors, a solution may result having no physical meaning whatsoever. One 
approach is to reduce the dimension of the problem so that A*A is computationally 
nonsingular. If a triangular decomposition of A 1 A is formed such that row operations 
are done in a selective order, this can be accomplished by examining the relative 
magnitudes of resulting diagonal elements as they are formed. This comparison is 
made relative to some a-priori significance level which allows for errors both in 
the data and in the computation. If the test fails, the column undtjr investigation is 
eliminated and the corresponding element of the solution vector X Q is assumed to 
be zero. 


The above procedure will result in a least squares solution to the linear system, 
but it has the disadvantage that observables are eliminated from the system, and full 
use is not made of the information available. The solution obtained is also rather 

arbitrary in that the elements of X fl which are assumed to be zero are determined by 
the order in which the reduction takes place. One way of circumventing these 
objections is to solve for the solution vector which has minimal Euclidian norm. 
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The problem can then be stated as follows: 


||AX 0 -B|| = min., 


|X 0 |) - min. 


The solution X^ exists and is unique regardless of whether A is of maximal rank. 


One way of obtaining the solution X Q is to make use of the Penrose pseudoinverse 
(see Appendix A), which has the following definition: For any rectangular matrix A, 
let A* denote the Penrose pseudoinverse of A. It can be shown that A* is well defined 
if it batisfies the following four axioms: 


The axioms imply that: 


(i) 

AA # A = A 

(ii) 

A # AA # - A # 

(iii) 

(A # A) T - A # A 

(iv) 

(AA # ) T = AA # . 

A # # 

= A, 

A # 

= A -1 if A” 1 exists 

<aV 

- (A ) . 


If we now consider the original linear system 


AX = B , 

and if we let 

X Q * A* B , 

then it follows that, for any vector X of dimension n , 

||A^-B|| < 1 1 AX - B 1 1 . 

Furthermore, if 

I|aV®H - H ax - b l f , 

then 

l|x 0 ||s l|x|| . 
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Hence, the Penrose pseudoinverse will provide the least squares solution with minimal 
Euclidean norm. 


LEAST SQUARES PROCEDURES INVESTIGATED 

Computational experiments were performed on the following five approaches to 
obtain a least squares solution. Three decimal digits of computational error were 
allowed for in the procedures that tested for reduction in computational rank. 

{1) A Gauss-Jordan matrix inversion algorithm applied to the normal equations 
(ref. 2 ). No check was made for computational singularity of A T A. The 

solution Xq was evaluated as follows: 

- T -1 T- 

X Q = (A 1 A) (A a B) . 

(2) An Andree pseudoinversion algorithm applied to the normal equations (see 

Appendixes B and C). Here a series of congruence transformations were made 
on the A^A matrix such that in each step the reduction was made to pivot about 
the largest remaining diagonal element. A computational reduction in rank was 
performed under the following test: 

Let be the first pivot element and c*. be a pivot element being tested. Then 
reduction in rank was performed when 



where k is the number of decimal figures carried in a computer word. Rank 
reduction was effected by setting and all subsequent pivot elements to zero. 

The solution X Q was evaluated as follows . 

X Q = (A a A f (A A B) . 


(3) A Gram-Schmidt orthoganalization procedure for evaluating the pseudoinverse, 
applied to the normal equations [see Appendixes B and E). Under this scheme 
an orthonormal basis for the column space of A^A was evaluated, and rank 
reduction was performed under the following test: 


Let 1 1 a j 1 1 represent the magnitude of a column vector of A 1 A being tested and 


let j | 1 1 be the magnitude of its projection on the independent column space 

already determined. Then reduction in rank was performed when 


115, - ?, II < 1(|3 


iM\ 

io k J 
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where k is the number of decimal figures carried in a computer word. The 
solution 5C was evaluated as follows: 


- T # T- 

X Q - (A 1 A) (A 1 B ) . 


(4) A Gram-Schmidt orthoganalization procedure for evaluating the pseudoinverse, 
applied to the rectangular m x n matrix A. The same routine as in (3) was used, 
but operations were done on the columns of A rather than A^A. The solution 
Xq was evaluated as follows: 


A # B 


(5) The Householder procedure for evaluating the least squares solution by finding 
the triangular decomposition of A T A directly from A (see Appendixes C and F). 
Under this scheme both A and B are pre-multiplied by an orthogonal matrix Q 
which transforms A into upper triangular form, with all elements below the 
main diagonal equal to zero; that is: 


«A - R - (f) • 



where R is an n x n upper triangular matrix and C 1 represents the first n 

elements of QB . Since the Euclidean norm is preserved under tn orthogonal 
transformation, 


| f AX - B || = } | RX - C 1 1 , 


and the least squares solution is clearly 


X 


0 



'V — 1 

provided that R exists. The triangulation of A was performed by operations 
on the columns of A in such an order that the diagonal elements formed were 
maximizrcl. Since R is the triangular decomposition of the positive definite 
matrix A 1 A, a test was made for computational reduction in rank by examining 
the relative magnitudes of the diagonal elements of R, Let r- be the first 
diagonal element formed and let r. represent a diagonal element being tested. 
Then R , and hence A , was assumed to be of reduced rank when 


r. < 10 3 

l 



* 


where k is the number of decimal figures carried in a computer word. 
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If rank was not maximal, no solution was evaluated. With maximal rank, X„ 
was evaluated by back substitution from the equation 



COMPUTATIONAL EXPERIMENTS 


In all of the experiments performed, attempts were made to recover the co- 
efficients for the polynomial Z 2 + 10 Z + 1 . Matrices were generated for situations 
of varying dimension and original data error. For a problem of dimension n, with 
m data points, the matrix A = (a..) and the vector B = (b.) were generated from 
a starting Z and A Z in the following manner: 1 


Z. = Z. , + AZ, i = 2, ••• , m 
i l-l 

b. = 1 + 10 Z. + Z.2 , i = 1 , ••• , m 

i ii’’ 

a.. = Z? 1 , j = 1, ••• , n 
ij i J 

i = 1, , m . 

— T 

If there were no computational errors involved, then the solution X„ = (1, 10, 1, 0, 0, • • •, 0) 
would clearly be the least squares solution, since 


|!ax 0 -b|| = o . 

Figures 1-a and 1-b represent the results with the matrices generated from 
= -1, AZ = 1/16, and m = 33 . The column dimension of A was varied from 5 to 

25. Single precision (24 binary digits) was used in the generation of matrices A and 
B; since each Z. can be exactly represented by four binary digits, no error was 

present in the A matrix up to column dimension 7. Results are presented with the 
least squares solution processes, including the formation of the normal equations,^ 
carried out both in single and double precision. In both instances matrices A and B 
were generated in single precision. Figure 1-a shows the common logarithm of the 

norm of the solution error vector (log.^ 1 1 X Q - (1, 10, 1, 0, • • • , 0) j | ) plotted 

against the column dimension of the A matrix. About 16 decimal digits are carried 
on the computer in double precision, and it is clear from tne graph that little sig- 
nificance was lost up to column dimension 7 for the two double precision procedures 
which avoid the formation of the normal equations. Figure 1-b represents the com- 
putational rank as determined by the various algorithms tested. The three routines 
which make use of the normal equations gave similar results in recovery until rank 
suppression occurred. In general, better recovery was observed with the Andree 
and Gram - Schmidt routines when rank was suppressed. Double precision 
arithmetic gave much better results up to dimension 20, at which point computational 
error became very high with the normal equations. It is evident that much better 
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results can be expected at high dimensions when the normal equations are avoided. 

Figures 2-a and 2-b demonstrate the results of data generated with high original 
error. Here matrices A and B were generated in single precision with row dimension 
m = 100, = 0.01 , and AZ = 0.01. Double precision was used to obtain the least 

squares"solutions. Most Z. generated here do not have exact binary representation, 
and a high degree of error was present in the higher powers of Z. . Up to dimension 
9, all procedures gave nearly the same results because the original data error in the 
A and B matrices was much greater than the computational error which occurred 
in the least squares soultions. Beyond dimension 9, the Andree algorithm applied 
to the normal equations gave the best recovery because there was greater suppression 
of computational rank. The Gram-Schmidt algorithm gave a different computational 
rank because non-negativeness was not forced in the (A T A)# matrix as it was with 
the Andree algorithm. Throughout the entire dimension range, both schemes which 
avoid the normal equations gave essentially the same results because original data 
error predominated. The results applying single precision to obtain the least squares 
solutions from these data matrices are shown in figures 3 -a and 3-b. At high, 
dimensions both single precision and double precision tended to give equally poor 
results, which further demonstrates the predominance of original data error. 

In figures 4-a and 4-b the same data matrices were used, but they were generated 
in double precision rather than in single precision. A marked improvement in re- 
covery is clearly evident, and the advantage of avoiding the normal equations is 
again demonstrated. The superiority of using the pseudoinverse with reduced com- 
putational rank when the normal equations are used is best shown in this example. 
Above dimension 11 the Andree algorithm gave the best recovery with the normal 
equations, and at very high dimensions results were comparable with the schemes 
which operate directly with the A matrix. However, at the high dimension level, 
computational noise became very great with all algorithms, so the comparison has 
little meaning at this point. 


COMPUTER RUNNING TIME 

Figures 5 and 6 show comparisons of computer running time in double and 
single precision for all schemes tested. The A matrix here is the one reflected 
in figures 2-a, 2-b, 3-a, and 3-b; that is, it was generated in single precision with 
m = 100, Z^ = 0.01, and AZ = 0.01. On the IBM 360/95 computer, which was used 
in all experiments, the minimum measurable time interval is 0.01664 seconds, so 
comparisons could not be made when the running time was less than this. Only the 
time required for the matrix inversion or pseudoinversion is reflected in these graphs, 
except for the Householder algorithm which does not evaluate an inverse. The back 
substitution time for obtaining the solution vector is included with the Householder 
scheme. 

In general, the Gauss-Jcrdan algorithm was the most economical in r unnin g 
time, and the Gram-Schmidt procedure applied to the A matrix required the most 
computer time. One can expect an even greater difference in running time if the 
A matrix has a larger row dimension. Observe that on the IBM 360/95 computer 
running time is essentially the same in both double precision and single precision 
wnere maximal rank is maintained. The Gram-Schmidt algorithm requires more 
time when there is reduction in rank, so a greater running time is shown for this 
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routine at the higher dimensions on the single precision plot. However, less time 
is required for the Andree algorithm with reduction in rank, so this routine appears 
more efficient in running time on the single precision graph at higher dimensions. 
The Householder algorithm also showed less running time with reduced rank because 
the solution was not completed. 


COMPUTER STORAGE REQUIREMENTS 

Storage restrictions may be a factor in choosing the subroutine one wishes to 
use in the least squares solution. All subroutines used in this study destroy the 
original data matrix, A or A T A, whichever applies. The table below lists the IBM 
360/95 storage requirements for all subroutines tested. The original matrix (A or 
A^A ) plus working arrays required are included in this list; with the the Householder 
algorithm the constant vector B is also included since it too is destroyed. The 
matrix A is assumed to be of dimension m x n. (Assume m = n for the Gram-Schmidt 
algorithm applied to the normal equations. ) 

T 

If m is large, one may be forced to store only the A A matrix and operate with 
the normal equations. Note that the Andree algorithm requires a large amount of 
working storage. 


Algorithm 

Double Precision 
Storage Requirements 
(Bytes) 

Single Precision 
Storage Requirements 
(Bytes) 

Gauss-Jordan 

9 Sn 

1876 + 8 (n + n + ^ ) 

1795 + 4 (n 2 4- 4n) 

Andree 

4538 + 8 (6n 2 + n) 

4472 + 4 (6n 2 + n) 

Gram-Schmidt 

2 

3304 + 8 (mn + n + 2n) 

3260 + 4 (mn + n 2 + 2n) 

Householder 

3190 + 8 (mn + m + 2n + ^ ) 

3150 + 4 (mn + m + 3n) 


CONCLUSIONS 

It is clear from these experiments that it best whenever possible to avoid 
formation of the normal equations in solving the linear least squares problem. 
However, storage restrictions may not allow the entire data matrix to be placed in 
core. In this case, some advantage is gained by making use of a pseudoinversion 
subroutine which allows for reduced computational rank. Of the two pseudoinversion 
schemes investigated, the Andree algorithm gave the best recovery when the normal 
equations were employed. This routine requires a great deal of temporary storage, 
however, and the Gram-Schmidt subroutine may be a necessary substitute. In some 
situations a Gauss-Jordan subroutine must be used if storage limitations are critical. 

Of the two procedures which operate directly on the original data matrix, the 
Householder algorithm is more efficient in both computer running time and storage 
requirements. However, it does not adequately handle situations of reduced com- 
putational rank. The Householder and the Gram-Schmidt routines gave comparable 
results in recovery. 
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There was a highly significant improvement in recovery when double precision 
rather than single precision was employed. This was true both with the least squares 
solution procedures and with the generation of the data matrices. Use of the normal 
equations with double precision generally gave better results than did the Gram- 
Schmidt and Householder schemes applied to the original rectangular data matrix in 
single precision. 
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COLUMN DIMENSION OF A-MATRIX 


z 

o 

55 

u 

w 

os 

ft 

w . 

ft! k 

cq w 

H 
ft) 

ft 

s 
o 
u 


8 

O 

z 

co 

ft) 


n ° 

0 C£> 

CO CO 

5 « 

< S 

1 a 

o 

o 

o S 

05 ^ 
PS 

w 

z 

o 

►— « 

H 
ft> 
ft! 

o 

co 


H 
ft) 
O 
05 
« 
ft 
CO 

z 
o 

co to 
W PS 
05 W 

II 

H 05 

S H 

a s 



12 


COLUMN DIMENSION OF A-MATRIX 
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COLUMN DIMENSION OF A-MATREX 
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Figure 3-o 


COMPUTATIONAL RANK COMPARISON FOR SINGLE PRECISION 
MATRIX INVERSION SUBROUTINES IBM 360/95 COMPUTER 
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COLUMN DIMENSION OF A- MATRIX 
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COLUMN DIi ' ENSiON OF A- MATRIX 
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COLUMN DIMENSION OF A -MATRIX 


RUNNING-TIME COMPARISON FOR F Gram-Schmidt 

DOUBLE PRECISION MATRIX INVERSION SUBROUTINES / Algorithm 

IBM 360/95 COMPUTER / (A-matrix) 



(*03S) 3WIX ONINNJ1H S6/09S WRI 


*o 

a> 

D 

U> 

\L 


18 


COLUMN DIMENSION OF A-MATRIX 
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COLUMN DIMENSION OF A MATRIX 
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APPENDIX A 


PENROSE PSEUDOINVERSE 


Let A # be defined as the Penrose pseudoimerse of A, where A # satisfies the 
four axioms: 


Axiom 1: 

AA # A 

= A. 


Axiom 2: 

A # AA # 

= A # . 


Axiom 3: 

(AA # ) t 

- A» t A t 

= AA # , 

Axiom 4: 

(A # A) T 

= a t a #t 

- A* A 


We will show that the Penrose pseudoinverse always exists and is unique. 

If the matrix A is non-singular, then the pseudoinverse is identical to the inverse. 

THEOREM A-l: If A -1 exists, then A # = A" 1 . 

Proof: 

AA*A = A, 

A' 1 AA # AA -1 - (A -1 A) A # (AA -1 ) - A -1 AA -1 , 

I A # I = A' 1 ! = I A" 1 = A- 1 , 

A* = A" 1 . 


THEOREM A -2: A # is unique. 


Proof: Assume X and Y are Penrose pseudoinverses of A. Then, both X 
and Y satisfy axioms 1 through 4. Thus: 


X -- XAX - (XA) X - (XA) t X = A T X T X - (A T Y T A T )X T X 


A-l 


= (A t Y t )A t X t X = (YA) t A t X t X = YAA t X t X = YA(A t X t )X 
= YA(XA) t X = YA(XAX) = YAX = Y(AX) = Y(AX) T = YX T A T 
= YX t A t Y t A t = Y(X t A t ) (Y t A t ) - Y(AX) t (AY) t 
= YAXAY - Y(AXA)Y - YAY = Y. 

THEOREM A-3: A #T = A T# . 

Proof: Since (A T ) # is the unique Penrose pseudoinverse of A T , it is sufficient 
to show that A #T satisfies axioms 1 through 4: 

(1) A A* A = A. 

Transposing, 

A t A* t A t = A t . 

(2) A # AA # = A # . 

Transposing, 

a» t a t a #t = a» t . 

(3) (AA») t = AA # = A #t A. 

Thus, 

(A ffT A T ) T = (AA # ) t = AA # = A #t A t . 

(4) (A # A) t = A t A #t = A # A . 

Thus, 

(A t A" t ) t = (A»A) t = A” A = A t A # t . 
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THEOREM A-4: (A # ) # - A. 


Proof: Since A* is the pseudoinverse of A, it follows that A is the pseudo- 
inverse of (A # ), from the symmetry of the axioms. By definition (A # ) # is the 
pseudoinverse of A # , thus (A # ) # = A by theorem A-2. 


THEOREM A-5: A* exists. 

Proof 1: Let A be a diagonal matrix: 

A = diag (\j, k r .... k n ) . 

Define A # to be the diagonal matrix: 

A# = diag {fx v /z 2 M n ) . 

where 

Mi = lA i , k. * 0 

Mi - 0 - 0 . 


We now show that A # satisfies axioms 1 and 2: 

(1) AA # A = diag (\. ) diag (/x. ) diag (A.) = diag^. 2 ^) = diag(\.) = A. 

(2) A # AA # = diag(^i) diag (/V. ) diag (/x. ) = diag^/V.) = diag(^) = A # . 


Since the product of diagonal matrices is diagonal, symmetry is preserved, thus 
satisfying axioms 3 and 4. 

Proof 2: Let A be symmetric; that is, let A = A T . Then, A has the repre- 
sentation 


A = S t DS, 
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where D is a diagonal matrix [D = diag (/V)] and S is orthogonal, that is S T = 
S -1 . Define A # as: 


A » = S t D # S. 


We now show that A # satisfies axioms 1 through 4: 


(1) 

AA # A 

- S T DSS T D# SS T DS# 

- S T DD*DS = S t DS = A. 

(2) 

A # AA # 

= S T D # SS T DSS t D # S 

= S t D # DD # S = S t D # S = A*. 

(3) 

> 

> 

3t 

H 

= (S t DSS t D # S) t = 

(S t DD*S) t = S t (DD # ) t S = s t dd*s 



= S T DSS T D*S - AA # . 

(4) 

(A#A) T 

= (S t D # SS t DS) t = 

(S t D # DS) t = S T (D # D) T S = S t D*DS 


- S t D # SS t DS = A n A . 

Proof 3: Let A be arbitrary, then we define A # as either 

A # = (A T A)# A t 
or 

A # = A t (AA t ) # , 

depending upon which resulting symmetric matrix, A t A or AA T , has the smaller 
dimension. 

We now show that A* satisfies axioms 1 through 4: 

(1) AA # A = A(A t A) # A t A 
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Let 


C ---• A(A t A) # A t A - A, 

then C T is given by 

C T = A t A(A t A) # A t - A T 

and C T C is given by 

C T C = |\ t A(A t A) # A t - A t J [\(A t A) # A t A - A 

- A t A(A t A) # A t A(A t A)* A t A - A T A(A T A)*A T A 
- A t A(A t A)* a t a + a t a 
= A T A - A T A - A T A + A t A 

= 0 . 

Thus, C = 0, or 

A(A t A)#A t A = A. 

(2) A*AA # = (A t A) # A t A(A t A)# A t = (A T A) # A T - A*. 

(3) (A*A) T = |(A t A)#A t aJ t = B *B j T = B # B = (A T A)#A T A = A # A. 

(4) (AA*) T = |a(A t A)* A t ] t = A t (A t A)# t A = A T (A T A) T# A 

= A T (A T A) # A = AA* 
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The main utility of the Penrose pseudoinverse is contained in the following result. 
Let the system of equations 


AX = B 


be given, and let 


— ♦ 

X 


o 


A*B. 


Then, for any X, 


ax 0 - b|| < Max - b|| . 


If for some X 


the equality holds, then 


liXo'll <11X11. 


Thus, the solution given by the Penrose pseudo inverse is the least squares solu- 
tion. If there is more than one least squares solution, the solution given by the 
Penrose pseudoinverse is the smallest solution in norm. 


THEOREM A-6 (Projection Lemma): 


|ax 0 - i! I 


MIN 

X 


I I AX - 5{ i 


if and only if 


X t A t (AX 0 - B) = 0 

for all X. 

Proof: Assume that 

X t A t (AX 0 - B) = 0 
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for all X , and let Y be given by 


Y = X 0 + Z. 

We then have 

| i AY “ B| ! 2 = j |AX 0 -B + AZ| t 2 - |!AX 0 -B|| 2 + 2 Z t A t (AX 0 -B)+ ||AZ|| 2 , 
since by assumption 

Z t A t (AXq - B) = 0. 

Then, 

1 |ay - b| I* = Max, -b|| 2 + ||az|| 2 . 

Thus, 

!|ay - B| | 2 < | jAX 0 - b| 1 2 

since 

llAZlI 2 >0. 

We now assume that 

| \A\ - B| ! < ItAY-Bil 

and wish to show that this implies that 

X T A T (AX 0 -B) = 0 

for all X, Assume that 


A-7 



X t A t (AX 0 -B) - a^O 
for some X. Let V and Y be defined as 


V 


-clX 

ilAXl! 2 


Y 


+ V. 


Then, 


May - b| i 2 - i I ax 0 - b| i 2 



AXa 

i I Axf | 2 


ilAXo - B| j 2 


- 2aX T A T (AX 0 - B) 

HaxT? 


I ! axI [ 2 

IlAXl! 4 


~ 2a 2 + a 2 _ ~a 2 

IlAXl! 2 | | AX| P " | | AX| 1 2 ’ 


which is a contradiction. 


THEOREM A-7: In the system of equations 


AX = B, 


where 

X 0 A«B 

and 


X i X 


o * 
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either 


I AX - B| | > | lAlJo - B| 


(AX - B| | = ||A^ -B|| and | |X| I > I IXqI! 


Proof: First we must show that 


X t A t (AX 0 - B) = 0 


as follows: 


X T A T (AA # B - B) = X T A T AA # B - X T A T B 


= X T A T (AA # ) B - X T A T B 


X T A t (AA^) t b - X T a t b 


= X T A T A #T A T B - X T A T B 


X T A T B - X T A T B 


Thus, from theorem A -6, 


l|AX-B|| > 1 1 AXjj - b| 


Assume the equality holds: 


Max -b|| - I|ax 0 -b! 
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If 


X 


x. + Y, 


then, by hypothesis, 


1 1 AX - Bl | 2 - I |AX„ - B| | 2 = 2Y t A t (AX 0 “ B) + ! |ay| | 2 


= llAYll 2 = 0. 


Now, 


llxjl 2 = I |x - (X - x 0 )| ! 2 = j |X| i 2 + | |x - x 0 [ I 2 - 2X t (X - X 0 ) 

= I |x| 1 2 + I |y|| 2 - 2 CX 0 + Y) T Y 

= ! |x| j 2 + | |y| j 2 - 2 B t a" t y- 2||y|| 2 

= I !x| I 2 - IlYll 2 - 2 B t A* t Y 

= I |x| 1 2 - I |y| I 2 - 2B t a* t a t a #t y 

= I |x| | 2 - I |Y| I 2 - 2 B t A* t (A # A) t Y 

= I |X| | 2 - llvll 2 - 2 B t A» t A t (AY). 

Since I I AY] | = 0, it follows that 

i |X 0 | | 2 = ||X|| 2 - ||Y|| 2 . 
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APPENDIX B 


PSEUDOINVERSE ALGORITHMS 


There exist many algorithms for the generation of the pseudoinverse of 
Penrose. Unfortunately, not much information is available in terms of their 
efficiency and computer requirements. Such a comparison is presently oeing 
made and will be the subject of a future report. In the interim, we will present 
two such algorithms which have been used extensively by the author. These 
algorithms possess the desirable property of allowing a limited control of the 
computational rank. 

In the first method, the pseudoinverse is formed by means of the Andree 
algorithm. This routine computes the pseudoinverse by means of the equation 

A # = (A T A) # A T 


Actually, the pseudo inverse of 


B = (A t A) 


is formed. This routine insures the symmetry of B # and also imposes the re- 
quirement that B # be non-negative. Thus, if computational errors leading to the 
loss of the positive definiteness of B exist, they are partially alleviated by this 
scheme. 

The second algorithm, which is basically a Gram-Schmidt orthogonalization 
procedure, can operate directly on the rectangular matrix A. There is no need 
to form the matrix associated with the normal equations, thus one major source 
of computational errors is eliminated. The penalty paid for this lies in the larger 
computer storage required to handle the rectangular matrix. 


ANDREE ALGORITHM 

The algorithm used in this paper Is a modification of the Andree algorithm 
by T. S. Englar 1 . The steps of this algorithm are ts follows: 

Kalman, R. D. and Englar, T. S., “An Automatic Synthesis Prog, ■'m for Optica' Filters and Control 
Systems,” NASA, July 1963. 
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1. Compute A T A or AA T , whichever has the smallest dimension. Call this 
resulting matrix B. It will be sufficient to compute B # , since A # is 
given by 


A* = (A T A) # A T or A # = A T (AA T ) # . 

2. Compute a non-singular matrix, S = (s L .), such that 

SBS t = E,’ 


where E = (e i .) is a diagonal matrix with elements either zero or one. 


3. If E = I then B is invertible, and 

B- 1 = S T S. 


If E ^ I then we define the matrix U = (lk .) by the following: 
For i ^ j, u. . - ~ s. . if e. . - 0, 

J i j i j ii ’ 


= 0 


For i - j , u . . - 0 

J ’ ii 

= 1 


if e. . 1. 

1 1 

if e. . - 0, 

if e. . - 1. 

1 1 


4. Compute 


C = U T BU . 


Delete the rows and columns of C corresponding to C . = 0, and call 
the resulting matrix D. Observe that D is a non-singular matrix of 
rank m. Compute D” 1 by means of the Andree algorithm as in step 2. 
Compute B # by means of 


B-2 


0 


D- 1 

B # - U ( 1 U T 

0 0 


5. Compute A # either by means of 


A# = B # A t 
or 

A # = A T B # . 


In the computation of step 2, the generation of the matrix S is done in 
at most 2n steps. The reduction is based upon pivoting in each step 
about the largest of the remaining diagonal elements. If after any step 
of iteration the largest of the remaining diagonal elements is less than 
the product of a preassigned constant, K, and the first pivotal element, 
then the rows and columns containing these elements are set to zero 
thus reducing the rank of the matrix. 


GRAM -SCHMIDT PROCEDURE 

This algorithm 2 permits one to pseudo invert a rectangular matrix, A, 
directly. The algorithm is based upon partitioning A in the form 


A = (R , RU), 


where all columns of R are linearly independent. The pseudoinverse A # is given 
by 


A* 


( (I + UU 7 )* 1 R # 

U T (I + UU 1 )" 1 R # 


2 Rust, B., Burros, W. R. and Schneeberger, C., ”A Simple Algorithm for computing the Generalized 
Inverse of a Matrix,” Communications of the ACM, Vol. 9, No. 5, May 1966, pp 381-387. 
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BS 




"U ■g|]Ji " ll !|P | 




This representation may be checked by substitution into the axioms. 

If (a p a 2 , ... a n ) is any set of linearly independent vectors, then it can be 
replaced by an orthonormal set (q v q 2 , ... c^) in the following manner: 



Observe that the resulting matrix Q is such that 

Q T Q = I- 


Let us apply this transformation to A = (R, RU) and keep track of the trans- 
formation by applying it simultaneously to the identity matrix partitioned as 


I 
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A will transform into 


A - (R , RU) ► (Q , 0) 


0 



Thus, we must have 


(R , RU) 


Z 

0 


\ = (RZ , RX + RU) = (Q , 0) 


or 


RZ = Q 


R = QZ" 1 


RX + RU = 0 


X - -U 


and R # is given by 


R # = Z Q t 

since 


R # = (R t R)- 1 R t = (Z- lT Q T QZ- 1 ) z- iT q t = ZQ t . 

Thus, the matrices R # and U are generated. It is still necessary to compute the 
matrices (I + UU T ) -1 and U T (I + UU T ) -1 . The first of these expressions can be 
written as 


(I + UU 1 )' 1 = I - U(U T U + I)-*U T . 
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If both sides are post -multiplied by (I + UU T ), the identity follows. The second 
expression can be written as 

U T (I + UU 1 )- 1 = (U T U + I)' 1 !! 1 , 

since 

U T (I + UU 1 )" 1 = U T - U T U(U T U + I)” 1 ^ 

= ^1 - U T U(U T U + I) -1 U T 

- |^I + (U T U + I)- 1 - (U T U + I)' 1 - U T U(U T U + I)- 1 U T 

= ^1 + (U T U + I)- 1 - (I + U T U) (I + U T 

= |l + (U T U + I)- 1 - I ju T 

= (U T U + 


If the Gram-Schmidt orthogonalization process is now applied to the matrix 

/“U \ 


this will transform to 
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where we have 


'-UPV /"UP 


P / \ P 


= P T U T UP + P T P = I , 


or this becomes 


Also, 


and 


P T (U T U + I)P = I 


(U T U + I) = P T_ 1 p - 1 


(U T U + I)' 1 = PP 1 


I - U(U T U + I)- 1 ^ = I-UPP T U T 


= I - (UP) (UP) T 


(U T U + I)- 1 ^ = P(UP) T 


Thus, takes the form 


A# - 


I - (UP) (UP) T j R # 
P(UP) T R# 
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APPENDIX C 


HOUSEHOLDER'S TRIANGULATION METHOD 
FOR OBTAINING THE LEAST SQUARES SOLUTION 


The following is a description of the triangulation scheme used in this study for 
obtaining the least squares solution directly from the data matrix without forming the 
normal equations. This procedure, originally described by Householder, has been 
developed into a workable algorithm for least squares problems by Golub. 1 The sub- 
routine and its description was taken from the IBM system 360 scientific subroutine 
package. 

In our linear system, 


AX - B , (1) 

both matrix A and vector B are premultiplied by an orthogonal matrix Q which trans- 
forms A into upper triangular form. All elements below the main diagonal will then be 
equal to zero; that is, 



where R is an n x n upper triangular matrix and C 1 represents the first n elements 
of QB . Since the Euclidean norm is preserved under an orthogonal transformation, 


| AX - B| | = 1 1 RX - C| | 


( 4 ) 


and the least squares solution is clearly 



provided S ^ exists. 


( 5 ) 


An effective way to realize this decomposition is via Householder transformations. 
The algorithm is a recursive n-step procedure defined by the following recursion 
formulas : 


A 


( 1 ) = 


V 1 ’) = A 


(6) 


^Golub, G. , "Numerical Methods for Solving Linear Least Squares Problems, " 
Numerische Mathematic, Vol. 7, No. 3 (1965), pp. 206-216. 

2 IBM System 360 Publication No. H20-0205-2, 1967, pp. 191-194. 
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A (k+1) . (k+1) ^(k) (k) , , o 

A' ’ - (a./ ) ' = P w A' k = 1, 2, • • • , n . (7) 

In order to get an upper triangular matrix A^ n+ every matrix P^ (k = 1, 2, • •• , n) 
should be defined so that it is symmetric and orthogonal, and so that it transforms all 
elements of the k th column of A'** below the main diagonal to zero. All of these 
restrictions to P^) are satisfied, setting 

P* k) - I -0 U (k) U (k) (8) 


^k / <k\ 

°k^k + a kk > 


™ (k) 2 I + for \k ' 0 

£ > <a J k) ) with 

i=k ) - for a^( k )< 0 


where I is the identity matrix. 


t j ' ’ denotes the transpose of column vector U v - (u. ') , the m components 


of which are defined as follows: 


u.' ’ =0 for i < k , 

i 


(k) , (k) 

v = 'k‘V - 


u. (k) = a.,^ for i>k. 

x lk 


Neither matrices P^ (k = 1, 2, • • • , n) nor matrix Q = P^ P^ n ^ • • • P^ are 
computed explicitly, since from (7) and (8): 


A (k + 1) - A (k) - U (k) 


Y <k) = ^ k ^ (k) A (k) . 


C-2 


— • /V\ 

Writing the components of row vector Y explicitly 


(k) 

V 

- 0 ; 

for j < k , 

(16) 

y (k) 

= 1 , 


(17) 

(k) 

y i 

^k 

f. forj>k . 

(18) 


Using (14), (16), (17), and (18), the explicit transformation formulas for 
matrix A' k ' appear as follows: 


\k 


<k + l) - 


- - ry 


k ’ 


(19) 


a. 


(k+1) 


ik 


= 0 , i = k + 1, k + 2, • •• , m , 


( 20 ) 


a.. (k + 1) = a.. (k) - u (k) y (k) , j = k+ 1, k + 2, • • • , n 

ij ij i j J ... 

i = k, k + 1, ••• , m 


( 21 ) 


(k ) 

All other elements of matrix A' ’ do not change. Naturaiiy the same transformation 
is performed on vector B . 


Using (6) and (7) 


B (1) - (b. (1) ) = B, 

(22) 

B^ + 1) - (b.^* 1 * = P (k) B (k) , k = 1,2, •••, n . 

(23) 

Using (£) and (23) : 

g(k + l) = g (k) j}(k) 

(24) 

with 

Z (k) = ^ k U (k) B (k) , 

(25) 
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or explicitly: 


with 


b ,( k + 1) = b.^-u.^ Z W 

1 J 1 


1,2, ••• , m 


(26) 


?(k) 






(27) 


In order to keep roundoff errors as small as possible, column interchange is 
performed in such a way that, at the stage, the column of A(^) is chosen to be 

reduced next which will maximize J a i c ] c ^ + ^| • Using (19) and (10), the index of 

this column is determined by giving the maximum overall j of: 


S j(k) = £ (a..^) 2 * J = K k + l> n 

i=k 


(28) 


After A^ + 1 ) has been computed, it is possible to compute s/ k + ^ as follows: 

s. (k + 1) = s. {k) - (a kj (k + 1) ) 2 , j = k + 1, k + 2, * * • , n , (29) 

since the orthogonal transformations leave the column lengths invariant. 

After having computed matrices A^ n+ ^ (upper triangular matrix) and B^ n + ^, 
the computation of the n by 1 solution vector X Q = (x.) is performed by back sub- 
stitution according to the following formulas: ! 


n ( n + 1) n 

nn 


b < n+1 > , 


(30) 


x,_ = 


*kk 


(n + 1) 


(v* 


1 ) 


- a 


(n + 1) 


(n + 1) 


k, k + 1 \ + l _a k, k + 2 *k + 2 


(n+1) 


(n+1) 


\,n-l x n- l~\n x 


\ k = n - j + 1 
n J ’ j = 2,3,..., n. 


(31) 
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After vector X. has been computed, back interchange of elements is performed 

according to column interchanges in matrices A v ' in order to get the correct 
sequence of components in the solution vector . 

The only case in which the whole procedure can fail occurs when, at any stage 
k, no column with nonzero parameter can be found; that is, no nonzero main 

diagonal element can be generated. In this case, the rank of matrix A is less than n. 
With respect to roundoff errors, the test is made in the following way. 

At first, all elements s^^ are computed according to (28), and the maximum 
overall j (i. e. , cj^) is determined. With the relative tolerance e given by input, the 
following absolute tolerance is generated: 


tol - ff, ( 


(32) 


If at any stage k (=1, 2, • • * , n + 1) the square root of the maximum overall 
j - k + 1, • ■ • , n of s ^ + ^ is not greater than cr » the rank of matrix A is 
declared to be k < n , and no solution can be found. 
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APPENDIX D 


ANDREE ALGORITHM SUBROUTINE 


The Andree algorithm for calculating the Penrose pseudoinverse of a matrix 
of equations has been converted into a FORTRAN subroutine. The listing for this 
subroutine appears on the following pages. 






SUBPOUTINE ANDPEECV.N.NR.EPS) 

DOUBLE PRFCI St-ON V.B.T * W t X * U »P . PR . DSQRT 

DIMENSION VI 30 , 30 ) ,Bi 30,30) ,R( 30) , T( 30.30) ,W< 30,30 ) ,X( 30.30 > 
DIMENSION 0(30,30) 

BILL = EPS 
NB = 56 
IND = 0 

LEE = 0 - ■ 

DO 333 I = i,N 

DO 333 J = i »N 

333 B( I , J) = V( I , J) 

66b DO 98 I = I , N 
P< I ) = 0 . 

DO 97 J s 1,N - 

97 rn.j) = o. 

98 TIM) i i, 

3b LEE = LEE+1 

IF CLEF-N) 77,77,78 
77 P = 0. 

~K = 0 

DO 22 I = 1 , N 

IF ( R ( I ) I 22.124,22 

124 IF (V(I,I) ) 82.82,83 

82 DO 84 J = t,N 

V ( I , J ) = 0. 

84 V(J,I) = 0. 

GO TO 22 

83 IF (P-V<T.l)) 21,22,22 

21 P = V ( I , I ) 

K = I 

IF (LEE.EQ.l) DM AX = P 

22 CONTINUE - - - - - - 

IF (P-I0.**BILL*DMAX/2.**NB) 28.23.7 

7 R< K> = K 

P = DSQRTC V( K ,K) ) 

DO 19 I = 1 ,N 
V< I . K) = V( I ,K)/P 

T< I » K ) = TH-.K1/P — - 

V(K, l ) = V ( K , I )/P 
19 T(K,l) = T(K,1)7P 
V( K, K ) = 1 . 

T(K.K) = I./P 
DO 25 I = I ,N 
IF <I-K) 26, 12b, 2b 

125 00 126 J = 1 ,N 

IF ( f-J) 127,128,127 

127 W< I , J) = 0. 

XC ! , J) « T< I , J) 

GO TO 126 

128 WT 1,1) » 1 • - • •• • 

X ( I , I ) = TCI. I) 

126 CONTINUE 
GO TO 25 

26 DO 1 0 J = 1 ,N 

«(I,J) = V( ! , J ) - V C I » K ) ♦ V ( K * J ) 

X(I,J) * T< I * 3 >— V < I ,K)*T (K » J > 
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10 CONTINUE 

a 5 continue 

DO '2* I = 1 , N 
DO 24 J = l.N 

V ( l « J ) = W ( I « J ) 

24 T< I, Jl * X< I i-J-» ■- - ■ - • 

GO TO 3b 

28 DO 30 ! =» l.N 

IF I »v C I > > 30.34.30 
34 DO 33 J- = l.N - - 

V( I . J) - 0. 

33 V< J, I > a 0. 

30 CONTINUE 

78 BILL = O. 

IF (IND) 38,39,38 

39 JOE = N 

DO 40 I = 1 ,N 
I r { V( f , I > ) 40.42,40 

42 JOE = JOE- l 

40 CONTINUE - 

IF (JOE-N) 43,44,43 

44 DO 45 I = l.N 
DO 45 II = l.N 
PR = 0. 

DO 46 J = l.N 

46 PR = RRtTU,II*TUiMF 

V ( I , I I > = PR 

45 CONTINUE 
GO TO 36 

43 DO 47 I = l.N 

IF ( R ( I ) ) 48.49,48 

49 DO 60 J « 1,-N — - 

60 U( I . J ) = — T ( I • J ) 

U< I « I ) = 0 . 

GO TO 47 

48 DO 61 J * l.N 

61 U< I . J) = 0. 

U( I , I ) •= 1~»- - - 

47 CONTINUE 

DO 50 I = l.N 
DO 50 II = l.N 
RR a 0. - - 

DO 51 J = l.N 

51 RR -= RR+Ul J , I-> +&{ J-il I 3 

W ( I , I 1 ) = RR 

50 CONTINUE - • 

DO 52 I = l.N 

OO 52 11 - l.N - 

RR = 0. 

DO 63 J » 4 ,N — *. - 

53 RR = RR + WC I . J) *U( J . I I ) 

Vdill) a RR 

52 CONTINUE 
IND = 1 
LEE - 0 
GO TO 666 



38 DO 34 ! = l.N 

DO 54 I! = t *N 
RR = 0 . 

DO 55 J = 1 »N 
55 RR = RR+T< J.I )*T( J.I I ) 
W{ 1.11} = PR 
54 CONTINUE 

HO 56 I = l.N 
DO 56 II = l.N 
RR = 0 . 

DO 57 J = l.N 

57 RR = RR+UI I .J)*W(J»I I ) 
5o T ( I . I I ) = RR 

DO 58 I = l.N 
DO 58 II = l.N 
RR = 0. 

DO 59 J = l.N 
59 RR = RR+T< I * J)*U< M vJ) 
V ( I . I I ) = RR 

58 CONTINUE 
36 NR = JOE 

RETURN 

END 



APPENDIX E 


GRAM-SCHMIDT PSEUDOINVERSION SUBROUTINE 


The Gram-Schmidt procedure for calculating the Penrose pseudoinverse of a 
matrix of equations has been converted into a FORTRAN subroutine. The listing 
for this subroutine appears on the following pages. 


SUBROUT INF GINV? t A,U, AFL A G . A TE MP , MR » NR , NC ♦ NR 1 .EPS) 
DOUBLE PRECISION A < MR , NC > . U < NC ♦ NC > . AFL AG t NC > , A TEMP< NC ) 
DOUBLE PRECISION F AC . DOT , DOT I , DO T2 • TOL . DSQRT 
DO 10 I = l.NC 
DO 5 J = 1 .NC 
5 O(I.J) = 0. 

10 U ( X • I ) = 1 . 

FAC = OOT< MR, NR, A. t . 1 > 

FAC = 1 • /DSQRT { F AC ) 

DO 15 I = l.NR 
15 A ( I , t ) = A ( I , 1 ) *FAC 
no 20 I = l.NC 
20 U ( I , 1 ) = U( I , 1 ) *FAC 
AFLAGtl) = l. 

N = 56 
NR I = NC 

TOL = ( 10,**EPS*.5**N> **2 

00 100 J = 2.NC 

DO T 1 = DOT < MR , NR , A , J , J ) 

JMl = J-l - 

DO 50 L = 1,2 

DO 00 K = l, JMl 

30 ATEMP(K) = DOT (MR*NR«A,J,K) 

DO 45 K = l, JMl 

DO 35 I = l.NR 

35 A ( I , J ) — M I. J)-ATEMP{K> *A< I ,K)*AFLAGtK> 

DO 40 I = l.NC 

40 U(I,J) = U< I . J)-AT£MPt K) *U< I .K> 

45 CONTINUE 
50 CONTINUE 

DO T 2 = DOTtMR.NR.A. J, J) 

IF { < OOT2/DGTI )-TOL) 5S-.55 t 70 
55 DO 60 I = 1 , JM 1 
ATEMPl I ) = 0. 

DO 60 K. = 1,1 

60 ATEMP(f) = A TEMP < 1 ) 4-U{ K. [ )*U<K,J) 

DO 65 l = l , NR 
A< I , J) = 0. - 

DO 65 K = l, JMl 

65 A ( I , J ) = At I . J)— A( I » K ) * A TEMPI K ) *AFLAG( K > 

AFLAG(J) = 0. 

FAC = DOT ( NC * NC , U , J , J ) 

FAC = 1 ./DSQRT(FAC) 

NR 1 = NR 1-1 
GO TO 75 
70 AFLAGtJ) = 1. 

FAC = I • /DSQR T t DOT 2 ) 

75 DO 80 I = l.NR 
00 At I , J ) = A( I , J ) *F AC 

DO 03 I * h NC - - - - 

05 Ut l, J) = U< I ,.J)*FAC 
100 CONTINUE 

DO 130 J = l.NC 
DO 130 I = l.NR 
FAC = 0. 

DO 120 K * J.NC 
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120 FAC ■= FAC+Ai I *K)*L'( J»K) 
130 A ( It.') = FAC 

RETURN 

CNO 
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DOUBLE PRECISION FUNCTION DOT ( MR ♦ NR . A . J . K > 
DOUBLE PRECISION A4MR,1),X 
X =' 0 • DO 
DO 50 I = 1 * NR 
X = X + A ( I ,J) **,( I « K ) 

50 CONTINUE 
DOT = X 
RETURN 
6 NO 



APPENDIX F 


HOUSEHOLDER’S ALGORITHM SUBROUTINE 


The Householder algorithm for obtaining a least squares solution is available 
as a FORTRAN subroutine in the IBM system 360 scientific subroutine package. The 
listing appears on the following pages. 


F-l 



n , . n a <*> o n n n <■*> r> 


\ 


SUBROUTINE OLLSQ (A»R*M*N*LtX. IPtV.EPS, IER.AUX) 

C — • - 

01 ME NS ION A< I ) ,B( 1 ) * X< t ) . IPIV< 1 ) , AUXl l ) 

DOUBLE PRECI SION A .+», X, AUX.P I V , H, S I G .BETA , TOL 
DOUBLE PRECISION DSORT 

C -ERROR- TE-ST- - 

IF ( M— N ) 30 » 1 « 1 

C GENERATION OF INITIAL VECTOR S K K 1.2.....N IN STORAGE 

C LOCATIONS AUX K K I«2,...,N 

1 PIV=0.D0 
lENO-O 
DO ♦ K= 1 . N 

IPIV<KI»K - - 

H=0.D0 

IST-IENO+I — ■- --- 

IEND= IENDTM 
OO 2 IMST.tENR 

2 H=H+A( I ) *A ( I ) 

AUXfK)»*4 

IF(H-PIV>4.4,3 

3 P I V*H 
KPl V=K 

♦ CONTINUE 

ERROR TEST - * - 

IF C P I V ) 31 * 3| ,5 

DEFINE TOLERANCE FOR CHECKING HANK OF A 
S SI G*DSQRT< PI V ) 

TOL=SIG*ABS(EPS) 


DECOMPOSITION LOOP- 

LM=-L* M 

ISTa-L - • 

DO 2* K= 1 »N 

I ST* rST'fKH - - - — - — - - • - - - - - 

1END= I ST+M-K 
I«KP I V-K 
IF ( I >8, 8, 6 

INTERCHANGE K-TH COLUMN OF A WITH KPIV-TH IN CASE KPIV.GT.K 

6 -H* AUXlKl - - - - - 

AUX(K)=AUX(KPIV) 

AUXIKPIV - — - 

ID=l*M 

DO 7 !«I ST. IENO 
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on nn non nn f*i n o n 


DO 10 1= ! ST « IEND 
10 - 
S 1 G= D SQR T ( S I G ) 

TEST ON SINGULARITY 

IF I S I G—TOL T32» 32 * -IT - - - — 

GENERATE CORRECT — SIGN-OF PARAMETER SI G 
1 1 H-=A C 1ST) 

IF 1+4 1 12* 1 3 « 13 

12 S I G— — S IG 

SAVE INTERCHANGE INFORMATION 

13 IPIV(KPIV»aIPfV<KT - 

IPIV(K)=KPl V 

GENERATION OF VFCTOR UK IN K-TH COLUMN OF MATRIX A AND OF 

PARAMETER SETA -- - ••• 

8ET'«-M rSIG 

A''SY> *OETA - •• - - -■ - - - 

RLTA=1 .D0/(SIG*BF rA) 

J»N+K _ _ 

AUX{ J ) = — SI G 

!F(K-N)U»r*H9 - 

TRANSFORMATION- YrF- M*TRTX A - - - 

14 PIV=0.D0 

- 10*0 ■ --- — - - — • 

JST=K+l 

KPf V* JST 

DO 18 J=JST,N 

- • — - 

H= 0 • DO 

DO 15 I* IS * •TEND 

I 1=1 AID 

15 H*H+A1 I ) *A ( I I I 

H=8FTA*H 

DO 16 1= 1ST c+ENO ~ ----- 

I I = I ♦ I O 

16 A< II >*A<II >-A<II*H 

UPDATING OF ELEMENT SJ STORED IN LOCATION AUX J - 
II=!ST*ID 

H=AUX< J>-AtTT F*AIHT — - 

AUXl J)=H 

iF<H-piv>ie,ie*»7 — - ----- 

17 PIV=H 

KPIV*J - - - - — -- - ■ ■ - - - 

18 CONTINUE 

C TRANSFORMATION OF RIGHT HAND SIDE MATRIX B 

19 DO 21 J*K.L»lrH 

H= 0 • D 0 

IENO=JAM-K — 

1 1 = I ST 

DO 20 Kh I€N» 
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on no 


H = H+ A { I i )*B( I ) 

20 1 1 *1 I ♦ I - - 

H=0ETA*H 

II*! ST 

DO 21 1 = J» ICND 



21 11 = 114-1 

C END OF DECOMPOSITION LOOP 

C 

C BACK SUBSTITUTION AND BACK INTERCHANGE 

I£R=0 - • • • - ■ 

I = N 

LN*L*N - - - 

PI V = 1 .D0/AUX( 2*N ) 

DO 22 K=N,LN,N ■ — - - - — * 

X< K)=PIV*BC I ) 

22 I=I+M 

IF (N- 1 ) 26* 26, 23 

23 JSTMN— 1 I+M+N - 

DO 25 J=2,N 

JST= JST-M— I 
K=N+N+1- J 

PI V= 1 .DO/AUMK ) 

KST=K— N 

ID=IPfVI KST >— -KfrT 

I ST=2— J 

DO 25 K= 1 rt_ — 

H=B( KST ) 

IST=I ST+N 
IEND= I ST+J-2 

I I - JST — 

DO 24 I = I S T * I END 
I I = I I CM - 

24 H=H— A ( II )*X { I ) 

1= f ST- 1 - 

11=1+10 

x( i >=xn i> 

XIII )=PI V*H 

25 KST=KST+M - - - 


COMPUTATION OF LEAST SQUARES 

26 !ST«N+1 — 

IEND=0 

DO 29 J»T*L 
I END= I END+M 
H=O.DO 

IF<M-N>29,29.27 

27 DO 28 I- 1ST, TEND - ■ 

28 H=HCB( I >*B( I ) 

I ST® I STEM 

29 AU X( J ) =H 
RETURN 

ERROR RETURN fN- CASE M LESS THAN N 
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on on 


30 IER=-2 
RETURN 

ERROR RETURN IN CASE OF ZERO-MATRIX A 

31 IE«=-1 

RETURN - - 

ERROR RETURN JN CASE— OF RANK OF MATRIX A CESS THAN N 

32 I ER = K— l 

RETURN - . 

END 



